14_blh_new-5km-shapefile

I ran the monthly maps code and noticed a couple differences, all are reasonable:

Not related to the change in shape files, but related to the map visualizations - We still need to decide if it makes more sense to facet based on fishing effort (where the VMS ping is, how it’s done now) vs. landings (where the ticket is landed, how it could be done). I think the answer should probably be similar to our decision to visualize by whale season, what makes the most sense from a management perspective? What regulations are the fishers held to, is it based on area fished or area landed?

I want to get a sense of the unique values and combinations of values for the new features, compared to the old features.

I’m going to load the two versions and compare.

grid_5km_v1 <- read_sf(here('spatial_data', 'master_5km_grid_tmer.shp'))
grid_5km_v2 <- read_sf(here('spatial_data', 'fivekm_grid_polys_shore_lamb.shp'))

What dimensions?

dim(grid_5km_v1)
[1] 76260    10
dim(grid_5km_v2)
[1] 72721    41

How many unique grid cells are there?

n_distinct(grid_5km_v1$GRID5KM_ID)
[1] 76260
n_distinct(grid_5km_v2$GRID5KM_ID)
[1] 72721

There are more grid cells in the older version, which makes me wonder if the extent is slightly broader. There are 31 more columns in the newer version, which makes sense because Blake made this to include more features I asked for.

What features do they each have?

colnames(grid_5km_v1)
 [1] "GRID5KM_ID" "centro_lat" "centro_lon" "STATE"      "WM_NGDC_M" 
 [6] "EEZ"        "WATER"      "Centroid_y" "Centroid_x" "geometry"  
colnames(grid_5km_v2)
 [1] "GRID5KM_ID" "ORIG_AREA"  "AREA"       "Shoreline"  "WM_NGDC_m" 
 [6] "NGDC_count" "NGDC_med_m" "NGDC_SD"    "NGDC_min_m" "NGDC_max_m"
[11] "NGDC_var"   "STATE"      "EEZ"        "Water_zone" "Salish_Sea"
[16] "RAMP_area"  "RAMP_zone"  "WDFW_WSMA"  "EEZ_MRGID"  "Derville"  
[21] "centro_lat" "centro_lon" "Centroid_y" "Centroid_x" "LL_EAST"   
[26] "UL_EAST"    "UR_EAST"    "LR_EAST"    "LL_NORTH"   "UL_NORTH"  
[31] "UR_NORTH"   "LR_NORTH"   "LL_LATITUD" "UL_LATITUD" "UR_LATITUD"
[36] "LR_LATITUD" "LL_LONGITU" "UL_LONGITU" "UR_LONGITU" "LR_LONGITU"
[41] "geometry"  
# what is the same?
intersect(colnames(grid_5km_v1), colnames(grid_5km_v2))
[1] "GRID5KM_ID" "centro_lat" "centro_lon" "STATE"      "EEZ"       
[6] "Centroid_y" "Centroid_x" "geometry"  
# what is unique to V1?
setdiff(colnames(grid_5km_v1), colnames(grid_5km_v2))
[1] "WM_NGDC_M" "WATER"    
# what is unique to V2?
setdiff(colnames(grid_5km_v2), colnames(grid_5km_v1))
 [1] "ORIG_AREA"  "AREA"       "Shoreline"  "WM_NGDC_m"  "NGDC_count"
 [6] "NGDC_med_m" "NGDC_SD"    "NGDC_min_m" "NGDC_max_m" "NGDC_var"  
[11] "Water_zone" "Salish_Sea" "RAMP_area"  "RAMP_zone"  "WDFW_WSMA" 
[16] "EEZ_MRGID"  "Derville"   "LL_EAST"    "UL_EAST"    "UR_EAST"   
[21] "LR_EAST"    "LL_NORTH"   "UL_NORTH"   "UR_NORTH"   "LR_NORTH"  
[26] "LL_LATITUD" "UL_LATITUD" "UR_LATITUD" "LR_LATITUD" "LL_LONGITU"
[31] "UL_LONGITU" "UR_LONGITU" "LR_LONGITU"

The newer version has more information about:

This is all based on the metadata file that Blake created, which is super helpful, since that detail is not knowable from column names alone.

What geometry type?

class(grid_5km_v1$geometry)
[1] "sfc_POLYGON" "sfc"        
class(grid_5km_v2$geometry)
[1] "sfc_MULTIPOLYGON" "sfc"             

Version 1 is polygons. Version 2 is multipolygons, because it trims the grid cells to the coast, and the shore can sometimes split a grid cell into multiple polygons.

What coordinate system?

st_crs(grid_5km_v1)
Coordinate Reference System:
  User input: WGS_1984_Transverse_Mercator 
  wkt:
PROJCRS["WGS_1984_Transverse_Mercator",
    BASEGEOGCRS["WGS 84",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ID["EPSG",6326]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["Degree",0.0174532925199433]]],
    CONVERSION["unnamed",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",31.96,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",-121.6,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",390000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]
st_crs(grid_5km_v2)
Coordinate Reference System:
  User input: CA_Curr_Lamb_Azi_Equal_Area 
  wkt:
PROJCRS["CA_Curr_Lamb_Azi_Equal_Area",
    BASEGEOGCRS["WGS 84",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ID["EPSG",6326]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["Degree",0.0174532925199433]]],
    CONVERSION["unnamed",
        METHOD["Lambert Azimuthal Equal Area",
            ID["EPSG",9820]],
        PARAMETER["Latitude of natural origin",30.5,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",-122.6,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["False easting",1000000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]

What extent?

st_bbox(grid_5km_v1)
     xmin      ymin      xmax      ymax 
-616809.7 -204655.9  979184.5 1929641.3 
st_bbox(grid_5km_v2)
   xmin    ymin    xmax    ymax 
  99280  -32131 1656582 2062869 

Well, they’re different, but that’s otherwise meaningless to me.

What’s it look like on a map?

background_map <- ne_states(country = 'United States of America')
bbox_v1 <- st_bbox(grid_5km_v1)
bbox_v2 <- st_bbox(grid_5km_v2)

ggplot() +
  geom_sf(data = grid_5km_v1, aes(fill = STATE), lwd = 0) + 
  geom_sf(data = background_map) +
  xlim(bbox_v1[1], bbox_v1[3]) +
  ylim(bbox_v1[2], bbox_v1[4])

ggplot() +
  geom_sf(data = grid_5km_v2, aes(fill = STATE), lwd = 0) + 
  geom_sf(data = background_map) +
  xlim(bbox_v2[1], bbox_v2[3]) +
  ylim(bbox_v2[2], bbox_v2[4])

What are the unique values for the new management zone variables?

unique(grid_5km_v2$RAMP_area)
[1] "Southern"      NA              "Central"       "Northern"     
[5] "San Francisco"
unique(grid_5km_v2$RAMP_zone)
[1] NA                "Zone 6"          "Zone 3 & Zone 7" "Zone 5"         
[5] "Zone 2"          "Zone 1"          "Zone 4"          "Zone 4 & Zone 7"
[9] "Zone 2 & Zone 7"
unique(grid_5km_v2$WDFW_WSMA)
 [1] NA      "57"    "60A-1" "61"    "58A"   "60A-2" "60D"   "59A-2" "58B"  
[10] "59A-1" "60C"   "59B"   "60B"  
unique(grid_5km_v2$Derville)
[1] NA                                       
[2] "North (North OR border to Cascade Head)"
[3] "South (Bandon to south OR border)"      
[4] "Central (Cascade Head to Bandon)"       

What do the EEZ & management zones look like on a map?

plot_grid <- function(grid, fill_var, background) {
  bbox <- st_bbox(grid)
  ggplot() +
    geom_sf(data = grid, aes(fill = get(fill_var)), lwd = 0) + 
    geom_sf(data = background) +
    xlim(bbox[1], bbox[3]) +
    ylim(bbox[2], bbox[4]) +
    labs(fill = fill_var)
}

plot_grid(grid = grid_5km_v2, fill_var = "EEZ", background = background_map)

plot_grid(grid = grid_5km_v2, fill_var = "RAMP_area", background = background_map)

plot_grid(grid = grid_5km_v2, fill_var = "RAMP_zone", background = background_map)

plot_grid(grid = grid_5km_v2, fill_var = "Derville", background = background_map)

plot_grid(grid = grid_5km_v2, fill_var = "WDFW_WSMA", background = background_map)

How many cells intersect the shoreline?

sum(grid_5km_v2$Shoreline)
[1] 1557

How many cells in the Salish sea?

sum(grid_5km_v2$Salish_Sea)
[1] 670

Are all the grid cells from the newer version, or are any lost?

length(setdiff(grid_5km_v2$GRID5KM_ID, grid_5km_v1$GRID5KM_ID))
[1] 0
length(setdiff(grid_5km_v1$GRID5KM_ID, grid_5km_v2$GRID5KM_ID))
[1] 3539

All of them are from the older version.

Where are the old grid cells?

old_grid_cell_ids <- setdiff(grid_5km_v1$GRID5KM_ID, grid_5km_v2$GRID5KM_ID)
old_grid_5km_v2 <- grid_5km_v1 %>% 
  mutate(is_old = GRID5KM_ID %in% old_grid_cell_ids) %>%
  filter(is_old)
old_bbox <- st_bbox(old_grid_5km_v2)

ggplot() +
  geom_sf(data = old_grid_5km_v2, color = 'orange') + 
  geom_sf(data = background_map) +
  xlim(old_bbox[1], old_bbox[3]) +
  ylim(old_bbox[2], old_bbox[4])

Seems like they’re all around the edges, which I won’t be using for these maps anyways, so no problems there.

I’m going to play around with some of the map aesthetics. How to rotate California coastline?

# non-rotated version
ca_grid <- grid_5km_v2 %>% filter(STATE == "CA")
plot_grid(ca_grid, "STATE", background_map)

plot(ca_grid, max.plot = 1)

# rotated version 1
# from https://r-spatial.github.io/sf/articles/sf3.html#affine-transformations
rot <- function(a) matrix(c(cos(a), sin(a), -sin(a), cos(a)), 2, 2)
ca_grid_rotated <- ca_grid
ca_grid_rotated$geometry <- ca_grid$geometry * rot(pi * 20/180)
plot(ca_grid_rotated, max.plot = 1)

# issue: there's no CRS when it's rotated this way, so geom_sf throws an error